# Bernd Beber, Philip Roessler, and Alexandra Scacco
# "Intergroup Violence and Political Attitudes: Evidence from a Dividing Sudan", Journal of Politics
# Code to generate figure A.3 in online appendix

# read data
    require(foreign)
    data <- read.dta("Beber_Roessler_Scacco_JOP_Data.dta")
    attach(data)
    sens <- read.dta("Beber_Roessler_Scacco_JOP_Sensitivity.dta")
    colnames(sens) <- gsub("_", ".", colnames(sens))
    attach(sens)

# plot to show simulation results from sensitivity analysis
    pdf("Beber_Roessler_Scacco_JOP_Sensitivity.pdf", width=8, height=8)
        # margins: bottom, left, top, right
        par(mar=c(5, 5, 2, 2))
        # distance to axis title, axis labels, axis line
        # par(mgp=c(2.4, .4, 0))
        plot(NA, xlim=c(-.4,.4), ylim=c(-.4,.4), xlab="Correlation with endogenous regressor", ylab="Correlation with outcome")
        # for red points, use `col=colors()[134]'
        points(.rho.d[.p>.05 | .coeff<=0], .rho.y[.p>.05 | .coeff<=0], pch=17)
        # for green points, use `col=colors()[104]'
        points(.rho.d[.p<=.05 & .coeff>0], .rho.y[.p<=.05 & .coeff>0], pch=1)
        # use `cluster.overplot' to minimize overlap?
        abline(h=0)
        abline(v=0)
        # add observed correlations
        # gender
        points(cor(female_r1, riots_affected_r1, use="pairwise.complete.obs"), cor(female_r1, separation_r1, use="pairwise.complete.obs"), pch="+", cex=1.4)
        # age
        points(cor(age_r1, riots_affected_r1, use="pairwise.complete.obs"), cor(age_r1, separation_r1, use="pairwise.complete.obs"), pch="+", cex=1.4)
        # working
        points(cor(working_r1, riots_affected_r1, use="pairwise.complete.obs"), cor(working_r1, separation_r1, use="pairwise.complete.obs"), pch="+", cex=1.4)
        # self-employed
        points(cor(selfemployed_r1, riots_affected_r1, use="pairwise.complete.obs"), cor(selfemployed_r1, separation_r1, use="pairwise.complete.obs"), pch="+", cex=1.4)
        # wealth index
        points(cor(wealth_factor_r1, riots_affected_r1, use="pairwise.complete.obs"), cor(wealth_factor_r1, separation_r1, use="pairwise.complete.obs"), pch="+", cex=1.4)
        # relative wealth
        points(cor(wealth_hh_hood_r1, riots_affected_r1, use="pairwise.complete.obs"), cor(wealth_hh_hood_r1, separation_r1, use="pairwise.complete.obs"), pch="+", cex=1.4)
        # education
        points(cor(educ_log_r1, riots_affected_r1, use="pairwise.complete.obs"), cor(educ_log_r1, separation_r1, use="pairwise.complete.obs"), pch="+", cex=1.4)
        # father's education
        points(cor(educ_father_log_r1, riots_affected_r1, use="pairwise.complete.obs"), cor(educ_father_log_r1, separation_r1, use="pairwise.complete.obs"), pch="+", cex=1.4)
        box()
    dev.off()

